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ABSTRACT 

We use a finite element approach based on Galerkin method to obtain approximate steady 
state solutions of the thermistor problem with temperature dependent electrical conductivity. 
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Q>^ \ 1 Introduction 

o 



In this paper we develop a method to approximate steady-state solutions of the following one- 
dimensional thermistor problem: 



O ■ ^"^ ^ ti ( \ \ ^t \\ \2 



dt dx (^(^)"^') = ^{u)\Vx\ , < j; < 1, t > 0, (1.1) 

^ ■ subject to boundary and initial conditions, 

■■■' k(u)— = -(3u ondnx(0,T), (1.2) 

ox 

u(x,0)=0, 0<x<l, (1.3) 

and coupled with the electric potential equation: 

{a{u)y^^)^ = 0, < X < 1, t > 0, (1.4) 

-^=^(x,t) onaO, (1.5) 

ox 

ip{x,0)=x, 0<x<l. (1.6) 
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The motivation for studying this kind of problem is that ()l.ip -( fL6]l has important imphcations 
for a variety of technological processes. For example, it arises in the analytical study of phe- 
nomena associated with the occurrence of shear band in metal being deformed at high strain 
rates [3j ; in the theory of gravitational equilibrium of polytropic stars [H] ; in the investigation 
of the fully turbulent behavior of flows [1]; in modelling aggregation of cells via interaction 
with a chemical substance (chemotaxis) [11]; and specially in modelling electrical heating in a 
conductor [H]. In this case, u is the temperature of the conductor, ip is the electrical potential. 
Functions a{u) and k{u) are, respectively, the electrical and thermal conductivities; P is the 
heat transfer coefficient. The condition (jl.3p is a condition of Robin-type. When /3 = it is 
called an adiabatic condition. Equation p.l|) consists in the heat equation with Joule heating 
as a source; (II. 4p describes conservation of current in the conductor. 

The thermistor problem has been extensively studied by several authors [H EJ El [7], where 
existence and uniqueness of solutions are given. Theoretical analysis, consisting in existence 
of solutions with the required regularity and which ensure error estimates of optimal order 
of convergence, are done in [8]. To construct a numerical approximation of the steady state 
solution we use a numerical method to approximate the solution of the parabolic problem. This 
approach has been used by [21 [10] in the one-dimensional thermistor problem. Further, in these 
last works authors consider the thermal conductivity k equal to 1 and a particular electrical 
conductivity, then they obtain the exact solution {(p{x, t) = x) of the conservation problem 
(|1.4|) - (|1.6|) and so system (|l.ip - (|1.6|) of thermistor problem is reduced to the following single 

heat conduction problem: 

du d'^u , , 

subject to the boundary conditions ()1.2p - (ll.3p . In this paper, we propose to solve both equations 
(jl.ip and (jl.6p at the same time by using a finite element method and a fully Crank-Nicolson 
approach. The formulation of the finite element method is standard and is based on a variational 
formulation of the continuous problem. In Section [2] we give the variational formulation of 
problem ()l.ip - ()1.6p . An algorithm for solving the problem is then proposed in Section [3l In 
Section JH numerical results are obtained for an appropriate test-problem. 

2 Variational formulation of the problem 

We divide the interval i7 = [0, 1] into A^ equal finite elements = xo<xi<...< xn = 1- Let 
(xj, Xj+i) be a partition of O and Xj+i — Xj = h = jj the step length. By S we denote a basis 
of the usual pyramid functions: 

\x+ (1 - j) on [x,j^i,Xj\, 

-\x+ {l+j) on [xj,Xj+i], 
otherwise. 

As indicated above, it is convenient to proceed in two steps with the derivation and analysis 
of the approximate solution of (|l.ip - (|1.6p . First, we write the problem in weak or variational 
form. We multiply the parabolic equation by Vj (for j fixed), integrate over (0,1), and apply 



Green's formula on the left-hand side, to obtain 

-^Vj dx + / k{u)VuVvjdx— I k{u)-—Vjds= / a{u)\V(p\ Vjdx. 

Using the boundary condition we get 

/ -—Vjdx+ / k{u)V uV V j dx + / (3uvj ds = / a{u)\Vip\ Vjdx. (2-1) 

Jn dt Jq Jqq Jq 

/ a{u)VipVvj dx = I a{u)-—Vjds. (2-2) 

Jo. Jdn (J^ 

We now turn our attention to solve this system by discretization with respect to the time 

variable. We introduce a time step r and time levels f„ = nr, where n is a nonnegative integer, 

and denote by vP' the approximation of u{tn) to be determined. We use the backward Euler 

Galerkin method, which is defined by replacing the time derivative in (|2.ip by a backward 

difference ^^^. So the approximations u""*"^, y?""'"^ admit unique representations 



We also have 



N N 



where a" , /x" are unknown real coefficients to be determined. Then, after decoupling, we 
have that 

/ ^^Vjdx+ / A;(n")Vu"+^Vz;j dx + / Pu''+^Vjds= / CTK)|Vvp'"|2i;j- dx, (2.3) 

Jn T Jn Jan Jn 



n -n+l 

aK)V(^"+iVi'jdx= / fiK)^-^^ ^jds. (2.4) 



and 

3 Formulation of the numerical method 

For scheme (12. 4p . we have 

AT 



1=1 "^^^ 






n+l n+1 

^ (a(u"(x,)) + a(n"(x,_i))) + ^ (a(n"(x,-+i)) + a(^."(x,_i))) 



n+l 

(a(n"(x,+i))+^(^"(x,))) 



/^^i 



2/i 



On the other hand, we have 



5^?"+^ f _ 

<y{u^) — T. — Vj ds = / a{u"')ipvjds 
an ov Jqq 

aK(l)MlH(l) -aK(0)MO)z;,(0) 

'-a(a^MO) ifj = 0, 
if j = l,...A^-2, 

ifj = iV-l. 

Using boundary conditions p.2p and initial condition (|1.3|) . it follows that 

W< ) J 



al^ = < 



Then, we have the resulting system of equations: 
for j = 0, 



n\\ , n+1 



(a(a^) + Ma-i) + 2^(ai )) f^o " (^(«-i) + 2^(«o) + ^(«i )) /"i 



-M0)(3a(a^)+a(a!!i)); (3.1) 



forj = l,...,iV-2, 



l^]ll {aia^) + a(a^_i)) + 2u«+i (a(a^+i) + a(a,"_i)) + ^"+^ (a(a^+i) + a(a,")) = ; 



(3.2) 



for j = N -I, 



= /i(^(l)(cT(a:^) + a(a^_i)). (3.3) 

Coming back to (j2.3p . the following may be stated in terms of the functions {vi)i: find the 
coefficients a^+^ in u"+i = J2f^-i <^7^^Vi such that 



V a^+^ / T;iT;j dx + r V a^+^ / A:(u")VuiVuj dx + t [ (iu^'^^Vj ds 
~^^ Jn ~^^ Jn Jon 

TV 

= V] a" / ViVjdx + T / o-(ii")|V(/?"pWj(ix. 



(3.4) 



In matrix notation, this may be expressed as 
where 



A = {aij) with element aij = / 



f jf j dx , 



B = (bij) with bij = / k{u^)VviVvj dx , 
Jn 



and 

a"'~^ is the vector of unknows (a^"^ )j=-i- 

Since the matrix A and B are Gram matrices, in particular they are positive definite and 
invertible. Thus, the above system of ordinary differential equations has obviously a unique 
solution. We solve the system ()3.ip for each time level. Estimating each term of (13. 1|) separately, 
we have: 

N N I 

y a^~^ / ViVj dx = 2, Q^r^ / ^i'^i '^^ 

= C(l-i / '^j-i'^j ^^ + (^1^ I ^j dx + a?+i / ^i^j+i d^ 

+ a'^ll / Vj-iVj dx + a"+^ / v] dx + v'j dx \ + a"+^^ / VjVj+i dx. 

Jxj^l \Jxj — l Jxj J Jxj 

Using the expression of Vj-i,Vj and Vj+i, we obtain 

X: «r^ / V.V, dx = ^^+1 + fa]^' + ^««:i. (3.5) 

In the same way, we have 

N „ N 



y <+! [ k{u'')VviVvi dx = y a"+^ [ k{u'')^^dx 
^, Jn .f^, Jn dx dx 



i=-l "" i=-l 



"?; £. ^-(""'^S "- + °r' £;■ M«")(^)= ^.- 



„+i r^+\,„.dvjdvj 



+1 



= ^/ k{u'')dx + ^-^ k{v^)dx-^-^ k{u!^)dx 

Q,"+l Q,"+l Q,"+l 

On other hand, we similarly have 

f a{u-)\^l\\dx=Y^ r^\{u^)WlW{x)dx 

i=i 
a«) 



It also holds: 



P 



f u^'+^j = /3'u"+i(l)fj(l) - /?u"+H0)vj(0) 

Jdn={o,i} 



-/3a^+i ifi = 0, 
if i = 1 . . . iV - 2, 

ifi = iV-l. 

Using together ()3.4p and ()3.5p , we get a system of A^ — 1 linear algebraic equations 

+ (^ - ^(^K+i) + ^(«i ))) «"+i - ^/3«r'^.(o) 



/i 



2/i 



/i 



n I ..n \2 



-a,"_i + -a] + -a^+i + -a{a]){-^,]_, + fi] + fi]+. 



Using the boundary conditions, we find 



a 



a_ 



n+l _ ^n+1 



a'r + 



a' + 



hf5 

hl3 



kK~') 



1 1 «r\ 



1 ) a^, 



,n+l 



A:(a^ 



N) 



a 



N 



/3h + A;(a^: 



,»i+l 



fe(a"-i^ 



/3h + A;(a^ 



]T"Af-l- 



From the initial condition we get 



Let 



ag = Oat = . 



^ - ^{kH) + Ho^-i)) ] , 



b=(Y + ^(fcK) + M«^)) 



and 



h T 



(M<) + A:(a^)) . 



6 2/i 

Substituting in (|3.6|) . we obtain the following system of equations: 
for j = 0, 



(3.6) 



1 + ^^^^) «o + ^«i + ^^K)(2a^o + mO)f ; (3.7) 



forj = l,...,iV-2, 



6~2h 






n+l 



2h 






for j = N -I, 






h 



h 



f ) «ri 

A; (a 



n-l^ 

AT 



6 6 \ /3/i + A;(a 



^ a^_i + -a(Q^_i)(2/x:^_i - ^^_2 + /i^(l))^ 



where 



AT 






/=(^-^(M«^) + M«^-i)))- 



4 An example 

In this section we give an example of a model of the thermistor problem: 

Ut = Uxx + l\'Px\'^ 

— — = — /3u on dQ 
ox 

^u{x,0)=0, 0<x<l 

= 
1 on dQ 



(4.1) 



(4.2) 



dip 
dx 
^ ip{x, 0) = X, < X < 1. 

The exact solution of the electrical potencial problem (j4.2p is ip{t, x)=x,0<x<l. Then, the 
diffusion equation (j4.ip can be reduced to the form 



du d'^x 



+ 7- 



dt dx^ 

Using the proposed Galerkin finite element approach, we get the following system of algebraic 

equations: 

for j = 0, 



h, 



Ph. 



h 



{aiiPh - 1) + &i - T/3) a^+^ + 2aia^+' = TT (1 + ^)«o + o < + l^h ; 



0.04 



0.03 



0.02 



0.01 




forj = l,...,iV 



Figure 1: The evolution of temperature. 



,n+l 



2h 



aia^+/ + bia''+' + aia^'jXl = -a^_i + — a^" + -a^+i + 7T/1 ; 



for j = N -I, 



aia^_2 +\h + 



oi 



/3/1+I 



a 



n+l 

N- 



h 



h 



l = 77«^-2 + 7T 4 + 



6 



6 



1 + /3/1 



ajv-i +lTh, 



where 



ai 



h T 2h 2t 



We now show some results from numerical experiments performed using our method and the 
computer algebra system Maple 10. According with physical situations, we choose values of (3 
and 7 verifying -g + 2 ^ -■ In particular, we fixed /3 = 0.2 and 7 = 0.1. The calculation of 
the steady-state for the thermistor problem is an important issue regarding the applications of 
the model in the thermistor device. We obtained stable steady-state times for r = 0.1, h = 0.01 
(see Fig. dj. 
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